//    Copyright (C) Kherty.  All rights reserved.
using System;
using OpenLS.Spreadsheet.Extensibility;

namespace OpenLS.Spreadsheet.AnalysisToolPak
{
    using Doub = Double;
    using Int = Int32;
    internal static class BesselFunctions
    {
        private static double BesselJ0(double x)
        {
            double j0;
            if (Math.Abs(x) < 8.0)
            {
                double y = x*x;
                double ans1 = 57568490574.0 +
                              y*
                              (-13362590354.0 + y*(651619640.7 + y*(-11214424.18 + y*(77392.33017 + y*(-184.9052456)))));
                double ans2 = 57568490411.0 +
                              y*(1029532985.0 + y*(9494680.718 + y*(59272.64853 + y*(267.8532712 + y*1.0))));
                j0 = ans1/ans2;
            }
            else
            {
                double ax = Math.Abs(x), z = 8.0/ax, y = x*x, xx = ax - 0.785398164;
                double ans1 = 1.0 +
                              y*(-0.1098628627e-2 + y*(0.2734510407e-4 + y*(-0.2073370639e-5 + y*0.2093887211e-6)));
                double ans2 = -0.1562499995e-1 +
                              y*(0.1430488765e-3 + y*(-0.6911147651e-5 + y*(0.7621095161e-6 - y*0.934945152e-7)));
                j0 = Math.Sqrt(0.636619772/ax)*(Math.Cos(xx)*ans1 - z*Math.Sin(xx)*ans2);
            }
            return j0;
        }

        private static double BesselY0(double x)
        {
            double y0;
            if (x < 8.0)
            {
                double y = x*x;
                double ans1 = -2957821389.0 +
                              y*(7062834065.0 + y*(-512359803.6 + y*(10879881.29 + y*(-86327.92757 + y*228.4622733))));
                double ans2 = 40076544269.0 +
                              y*(745249964.8 + y*(7189466.438 + y*(47447.26470 + y*(226.1030244 + y*1.0))));
                y0 = (ans1/ans2) + 0.636619772*BesselJ0(x)*Math.Log(x);
            }
            else
            {
                double z = 8.0/x, y = z*z, xx = x - 0.785398164;
                double ans1 = 1.0 +
                              y*(-0.1098628627e-2 + y*(0.2734510407e-4 + y*(-0.2073370639e-5 + y*0.2093887211e-6)));
                double ans2 = -0.1562499995e-1 +
                              y*(0.1430488765e-3 + y*(-0.6911147651e-5 + y*(0.7621095161e-6 + y*(-0.934945152e-7))));
                y0 = Math.Sqrt(0.63661972/x)*Math.Sin(xx)*ans1 + z*Math.Cos(xx)*ans2;
            }
            return y0;
        }

        private static double BesselJ1(double x)
        {
            double j1;
            if (Math.Abs(x) < 8.0)
            {
                double y = x*x;
                double ans1 = x*
                              (72362614232.0 +
                               y*
                               (-7895059235.0 + y*(242396853.1 + y*(-2972611.439 + y*(15704.48260 + y*(-30.16036606))))));
                double ans2 = 144725228442.0 +
                              y*(2300535178.0 + y*(18583304.74 + y*(99447.43394 + y*(376.9991397 + y*1.0))));
                j1 = ans1/ans2;
            }
            else
            {
                double ax = Math.Abs(x), z = 8.0/ax, y = z*z, xx = ax - 2.356194491;
                double ans1 = 1.0 + y*(0.183105e-2 + y*(-0.3516396496e-4 + y*(0.2457520174e-5 + y*(-0.240337019e-6))));
                double ans2 = 0.04687499995 +
                              y*(-0.2002690873e-3 + y*(0.8449199096e-5 + y*(-0.88228987e-6 + y*0.105787412e-6)));
                double ans = Math.Sqrt(0.636619772/ax)*(Math.Cos(xx)*ans1 - z*Math.Sin(xx)*ans2);
                j1 = (x < 0.0) ? -ans : ans;
            }
            return j1;
        }

        private static double BesselY1(double x)
        {
            double result;
            if (x < 8.0)
            {
                double y = x*x;
                double ans1 = x*
                              (-0.4900604943e13 +
                               y*
                               (0.1275274390e13 +
                                y*(-0.5153438139e11 + y*(0.7349264551e9 + y*(-0.4237922726e7 + y*(0.8511937935e4))))));
                double ans2 = 0.2499580570e14 +
                              y*
                              (0.4244419664e12 +
                               y*
                               (0.3733650367e10 + y*(0.2245904002e8 + y*(0.1020426050e6 + y*(0.3549632885e3 + y*1.0)))));
                result = (ans1/ans2) + 0.636619772*(BesselJ1(x)*Math.Log(x) - 1.0/x);
            }
            else
            {
                double z = 8.0/x, y = x*x, xx = x - 2.356194491;
                double ans1 = 1.0 + y*(0.183105e-2 + y*(-0.3516396496e-4 + y*(0.2457520174e-5 + y*(-0.240337019e-6))));
                double ans2 = 0.04687499995 +
                              y*(-0.2002690873e-3 + y*(0.8449199096e-5 + y*(-0.88228987e-6 + y*0.105787412e-6)));
                result = Math.Sqrt(0.636619772/x)*(Math.Sin(xx)*ans1 + z*Math.Cos(xx)*ans2);
            }
            return result;
        }

        public static double BesselY(double x, int n)
        {
            if (n == 0) return BesselY0(x);
            if (n == 1) return BesselY1(x);
            if (x <= 0) throw new ArgumentException("x");
            double tox = 2.0/x, result = BesselY1(x), bym = BesselY0(x);
            for (int i = 1; i < n; ++i)
            {
                double byp = i*tox*result - bym;
                bym = result;
                result = byp;
            }
            return result;
        }

        public static double BesselI(double x, double order)
        {
            if (x < 0)
            {
                if (order != Math.Floor(order))
                    throw new ErrorException(Errors.NumError);
                if (order%2 == 0)
                    return BesselI(-x, order);
                return -BesselI(-x, order);
            }
            if (order == 0)
                return BesselI0(x);
            if (order == 1)
                return BesselI1(x);
            throw new NotImplementedException();
        }

        private static double BesselI0(double x)
        {
            double positiveX = Math.Abs(x);
            if (positiveX < 3.75)
            {
                double y = x/3.75;
                y *= y;
                return 1.0 +
                       y*(3.5156229 + y*(3.0899424 + y*(1.2067492 + y*(0.2659732 + y*(0.360768e-1 + y*0.45813e-2)))));
            }
            else
            {
                double y = 3.75/positiveX;
                return (Math.Exp(positiveX)/Math.Sqrt(positiveX))*
                       (0.39894228 +
                        y*
                        (0.1328592e-1 +
                         y*
                         (0.225319e-2 +
                          y*
                          (-0.157565e-2 +
                           y*(0.916281e-2 + y*(-0.2057706e-1 + y*(0.2635537e-1 + y*(-0.1647633e-1 + y*0.392377e-2))))))));
            }
        }

        private static double BesselI1(double x)
        {
            double positiveX = Math.Abs(x);
            if (positiveX < 3.75)
            {
                double y = x/3.75;
                y *= y;
                return x*
                       (0.5 +
                        y*
                        (0.87890494 +
                         y*(0.51498869 + y*(0.15084934 + y*(0.2658733e-1 + y*(0.301532e-2 + y*0.32411e-3))))));
            }
            else
            {
                double y = 3.75/positiveX;
                return
                    Math.Abs((Math.Exp(positiveX)/Math.Sqrt(positiveX))*
                             (0.39894228 +
                              y*
                              (-0.3988024e-1 +
                               y*
                               (-0.362018e-2 +
                                y*
                                (0.163801e-2 +
                                 y*
                                 (-0.1031555e-1 +
                                  y*(0.2282967e-1 + y*(-0.2895312e-1 + y*(0.1787654e-1 - y*0.420059e-2)))))))));
            }
        }


        public static double BesselJ(double x, int n)
        {
            if (x == 0.0) return 0;
            if (n == 0) return BesselJ0(x);
            if (n == 1) return BesselJ1(x);
            double rj;
            const int iacc = 40;
            const double bigno = 1.0e10;
            const double bigni = 1.0e-10;
            if (Math.Abs(x) > n)
            {
                double tox = 2.0/Math.Abs(x), bjm = BesselJ0(Math.Abs(x)), bj = BesselJ1(Math.Abs(x));
                for (int j = 1; j < n; ++j)
                {
                    double bjp = j*tox*bj - bjm;
                    bjm = bj;
                    bj = bjp;
                }
                rj = bj;
            }
            else
            {
                double tox = 2.0/Math.Abs(x);
                int m = 2*((n*(int) Math.Floor(Math.Sqrt(1.0*(iacc*n))))/2);
                double ans = 0, jsum = 0, sum = 0, bjp = 0, bj = 1.0;
                for (int j = m; j > 0; --j)
                {
                    double bjm = j*tox*bj - bjp;
                    bjp = bj;
                    bj = bjm;
                    if (Math.Abs(bj) > bigno)
                    {
                        bj *= bigni;
                        bjp *= bigni;
                        ans *= bigni;
                        sum *= bigni;
                    }
                    if (jsum != 0) sum += bj;
                    jsum = 1 - jsum;
                    if (j == n) ans = bjp;
                }
                sum = 2*sum - bj;
                rj = ans/sum;
            }
            if ((x < 0.0) && ((n%2) == 1)) rj = -rj;
            return rj;
        }

        internal static object BesselK(double x, int n)
        {
            if (n < 0)
                throw new ErrorException(Errors.NumError);

            if (x == 0) return BesselK0(x);
            if (x == 1) return BesselK1(x);

            double t = BesselK0(x);
            double result = BesselK1(x);
            for (int i = 1; i < n; i++)
            {
                double temp = t + i*2/x*result;
                t = result;
                result = temp;
            }
            return result;
        }
      
        private static double BesselK1(double x)
        {
            if (x <= 2.0)
            {
                double x2 = x*0.5;
                double y = x2*x2;
                return Math.Log(x2)*BesselI(x, 1) +
                       (1.0 +
                        y*
                        (0.15443144 +
                         y*(-0.67278579 + y*(-0.18156897 + y*(-0.1919402e-1 + y*(-0.110404e-2 + y*(-0.4686e-4)))))))/x;
            }
            else
            {
                double y = 2.0/x;
                return Math.Exp(-1*x)/Math.Sqrt(x)*
                       (1.25331414 +
                        y*
                        (0.23498619 +
                         y*(-0.3655620e-1 + y*(0.1504268e-1 + y*(-0.780353e-2 + y*(0.325614e-2 + y*(-0.68245e-3)))))));
            }
        }

        static double bessk0(double x)
        /*------------------------------------------------------------*/
        /* PURPOSE: Evaluate modified Bessel function Kn(x) and n=0.  */
        /*------------------------------------------------------------*/
        {
            double y, ans;

            if (x <= 2.0)
            {
                y = x * x / 4.0;
                ans = (-Math.Log(x / 2.0) * BesselI0(x)) + (-0.57721566 + y * (0.42278420
                   + y * (0.23069756 + y * (0.3488590e-1 + y * (0.262698e-2
                   + y * (0.10750e-3 + y * 0.74e-5))))));
            }
            else
            {
                y = 2.0 / x;
                ans = (Math.Exp(-x) / Math.Sqrt(x)) * (1.25331414 + y * (-0.7832358e-1
                   + y * (0.2189568e-1 + y * (-0.1062446e-1 + y * (0.587872e-2
                   + y * (-0.251540e-2 + y * 0.53208e-3))))));
            }
            return ans;
        }

        private static double BesselK0(double x)
        {
            return bessk0(x);
            if (x <= 2)
            {
                double x2 = x*0.5;
                double y = x2*x2;
                return -1*Math.Log(x2)*BesselI(x, 0)
                       +
                       (-0.57721566 +
                        y*
                        (0.42278420 + y*(0.23069756 + y*(0.3488590e-1 + y*(0.262698e-2 + y*(0.10750e-3 + y*0.74e-5))))));
            }
            else
            {
                double y = 2.0/x;
                return Math.Exp(-1*x)/Math.Sqrt(x)*
                       (1.25331414 + y*(-0.7832358e-1 + y*(0.2189568e-1 + y*(-0.1062446e-1 + y*
                                                                             (0.587872e-2 +
                                                                              y*(-0.251540e-2 + y*0.53208e-3))))));
            }
        }
    }
}